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Abstract 

A simple and efficient method for quantum Monte Carlo simulation is presented, based 
on discretization of the action in the path integral, and a Gaussian averaging of the 
potential, which works well e.g. with the Coulomb potential. Non-integrable hard core 
potentials can not be averaged in the same way. 
PACS numbers: 02.70.-c, 02.70.Ss, 05.10.-a, 05.30.-d 
Keywords: quantum Monte Carlo, partial averaging, Coulomb potential. 

1 Introduction 

The quantum Monte Carlo method is well established as an efficient calculational tool for 
many-body problems. See for example the review articles It is well suited for bosonic 

systems without magnetic field, where the path integral has only positive contributions. But 
also femionic systems can be treated, in spite of the troublesome "sign problem", just one 
application is the computation of the high temperature phase diagram of hydrogen j3J Oj]. 
One limitation is that it is a statistical method, so that every factor of 10 in precision costs 
a factor of 100 in computing time. But it does not hit the "exponential wall", because it 
only needs to represent particle positions, or particle paths, and so the number of parameters 
increases linearly with the number of particles. 

By contrast, the number of particles that can be handled by methods based on computing 
realistic manybody wave functions is limited by the exponential increase in the number of 
parameters needed for describing such wave functions. The density functional method [0] is 
less severely limited, since it uses one particle wave functions, but on the other hand it has 
to rely on clever approximation techniques. 

The purpose of the work presented here was to look for a simple and efficient way of 
handling the Coulomb potential in quantum Monte Carlo simulation. Formulae for the exact 
propagator of the two-particle Coulomb problem are known, and can even be derived by 
path integral methods IE] > but it is not clear whether they are useful for the simulation of 
many-body systems. In ref. [3] fitted formulae for two-particle propagators were used. This 
may be a good enough method, but it may nevertheless be of interest to look for more direct 
approaches. 

The Fourier representation of paths in path integrals was introduced by Feynman, to- 
gether with the idea of approximately integrating over infinitely many Fourier components 
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by averaging the potential 9 . The method was further developed by Doll, Coalson, and 
Freeman |1U1 lllj under the name of partial averaging. In their work and in subsequent work, 
see e.g. ^21 E]) the object focused upon has been the propagator, involving paths from one 
point to another, more than the partition function, involving closed paths. 

The basic idea is to let the path integral include only paths represented by finite Fourier 
series with a fixed number of terms, and to think of each such path as representing all the 
infinite Fourier series to which it can be extended. The result is that, as long as correlations 
along the path are neglected, each point on one truncated path represents a Gaussian dis- 
tribution of points, with a standard deviation decreasing from the middle of the path and 
vanishing at the end points. In the case of the Coulomb potential, the Gaussian averaging has 
the important effect of removing the singularity at zero distance. It is then a complication 
that the averaging varies along the path, disappearing towards the end points |13j . 

The modification proposed here is to average with a standard deviation which is constant 
along the path. This seems a natural approach when the partition function is computed 
directly, and not via the propagator. In more detail, the method proposed amounts to a 
discretization of the action integral, with an averaged potential, and the computation of the 
kinetic energy part of the action by means of a finite Fourier transform. As discussed below, 
within this method it is easy to add to any potential an auxiliary confining harmonic oscillator 
potential, which makes the partition function mathematically well defined. 

A different topic which is not addressed here is the optimization of the Monte Carlo 
sampling procedure. See in this connection the comment in an appendix of ref. ^3]. A finite 
Fourier transform is a central part of the present method. In the method as formulated here 
it is assumed that the number of time steps is odd, hence the standard fast Fourier transform 
with 2 n points, n = 1, 2, . . ., can not be used. One solution is to use 3 n points and the fast 
Fourier transform to base three However, it would also be straightforward to modify the 
method so as to use an even number of time steps. 



2 The imaginary time path integral 

Equilibrium properties of a physical system at a finite temperature T can be computed from 
the partition function 

Z(J3) = Tre-P H , (1) 

where H is the Hamiltonian, and (3 = 1/(A;^T). One may regard j3H formally as an imaginary 
time interval. 

To be specific, we consider most of the time one particle of mass m in three dimensions. 
The Hamiltonian is H = T + V, with T = p 2 /(2m) the kinetic energy and V = V(r) the 
potential energy. The partition function has the following path integral representation, 



Z = C 3 [ d 3 a J] fe? / d 6 a n \ exp 

n=l 



(2) 



We define the constants 



„ 1 / m „ 2nir frn „ . , 
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and C n = Cnj^pn for n = 0, 1, 2, 5 is the imaginary time action, 

S 1 r?* , /l o \ ^ „o, ., l 



ft ft 



^ dr^-m(Hr)) 2 + V(r(T))j=Y,C n 2 \a n \ 2 + -j o dr 



V(r(r)) . (4) 



The path r = r(r) is periodic in the imaginary time r, with period /3ft, and is given by the 
infinite Fourier series 



r(r) = 2, a»e' . (5) 

n=— oo 



The Fourier components a n are complex and satisfy the relations a_ n = a^, so that r(r) is 
real. In particular, ao is real. The time derivative r is with respect to the imaginary time r. 

This Fourier expansion is the natural one in the computation of the partition function, 
which involves periodic paths. A slightly different expansion is needed in the computation of 
propagators, see e.g. [31 ITU1 ITTl IT21 ITS] . 

It may be useful to sketch the derivation of this Fourier path integral. We start from the 
approximation 

Z^Trf[(e-i T e-7 y ). (6) 

3=1 

We insert J times the identity operator / = Jd 3 r |r)(r| where \r) is the position eigenstate, 
and introduce the free particle propagator 

(r'\e-f 3T \r) = C 3 e- c 'o\ r -r'\\ (7) 

to obtain the approximation 

Z« (7jQ,) 3J / d 3 n d 3 r 2 • • • d 3 r j exp(-^\ . (8) 
Here Sp is the "primitive" discretized action defined by 

3=1 3=1 

We define rj+\ = r\. This approximate expression for Z is exact when V = 0, and gives 
then in particular for J = 1 that 

Z = C 3 Jd 3 r 1 . (10) 

To make this integral finite we should regularize, e.g. by introducing periodic boundary con- 
ditions or an external harmonic oscillator potential. See Section below. 

We now take J to be odd, J = 2K + 1, and make the finite Fourier transform 

K „ . 

rj = }2 fl « e 2K+1 ■ ( n ) 

n=-K 
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It gives that 



q oo a J 

^ = E^ 2 xKI 2 + 5£^)> (12) 

71 = 1 j = l 



with, for n = 1,2,..., 



2(2if + 1) fm . ( rm \ ~ C n>K 



In order to transform the integral over the positions r,- into an integral over the Fourier 
components a n , we note that 



£ |dr/ = (2K + 1) (\da \ 2 + 2 £ |da n | 2 ] , 

7 V n=l / 



(14) 



and hence, 

A' 



J] d 3 r, = 2 3K (v / 2KTT) 3(2 ^ +1) d 3 a ]J d 6 a n . (15) 

j 71=1 



Using the identity 

sin I 

2if + 1 



n=l 

we rewrite Equation (jSJ as follows, 

A" 



Z « C 3 1 d 3 a f f[ C n 6 ^ | d 6 a n j exp (-^) . (17) 

The limit K — > oo gives Equation (J2J). We regard as a function of the variable r = 
j(3fi/{2K + 1), which becomes continuous in the limit. 



3 Averaging the potential 

The integrand of the path integral, Equation (j2j), is the negative exponential of S/h, Equa- 
tion Q). One way to interpret this is that the kinetic part of S defines independent Gaussian 
probability distributions for the Fourier coefficients a n , n > 0, such that the real and imagi- 
nary parts of the x, y, z components of a n have mean values zero and standard deviations 

1 h 
an ~ V2C n ~ 2^ 

The path integral may be computed approximately by the partial averaging method |10[ 
lllj . We integrate explicitly over the lowest Fourier components, and approximate the integral 




4 



over the infinite number of remaining coefficients simply by averaging the potential. This 
means that we choose some finite K and define 



K 

2riTTT 



R(t) = a n^~ ■ (19) 



n=-K 

The remainder term s(r) = r(r) — R(t) will have a Gaussian distribution with zero mean 
and with variances 

(Mr)) 2 ) = ((s y (r)) 2 ) = (Mr)) 2 ) = = a 2 . (20) 

To compute a, we compute 

(s(n) ■ s{t 2 )) _ (in 2 ~ 1_ / 2mr(n-72) \ _ /3ft 2 

3 " 2vr 2 m n= ^ +1 n2 COS v /ffi J ~ 2vr 2 m /Jft J' 1 j 

where the function fx = fx{u) has period 1 in its argument u, 

„ , , ^ cos(2n7rii) 9 / 1\ 2 7r 2 cos(2n7rn) 

^ 1,-^=" (—2) ~ 12 ~ ^ n 2 ' 

n=K+l n=l 

The last formula is valid for < u < 1. Thus we have 



(22) 



~ 2vr 2 m /k(0) 2vr 2 m \ 6 ra 2 J ~ (2# + l)^ 2 m ' (23) 

introducing the approximation 

00 1 f°° dn 2 

which is about 20% larger than the exact result 7r 2 /6 in the worst case K = 0. 
The problem facing us is to compute the integral 

1= [pnJ^o, n exp(-C n 2 |a„| 2 )) exp^-i^drF(r(r)^ . (25) 

Here r(r) = -R(t) + s(r) is given by Equation (JSJ), and R(r) by Equation (|T§|) . To simplify 
our notation, we write the integral as an average. Next, we assume that s(t\) and s(t2) are 
uncorrelated for n 7^ T2, which is true to a certain approximation, as shown by Equation (|21jl. 
In this approximation we may compute the integral by the following formal reasoning, 

/ = (exp^-y^\rV(r(r))\\ = / jj eX p(-^ y(r(r))) \ 

« T V(r(r)))=exp -- ^ drW(R(r))j, (26) 
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where W is an averaged version of the potential V, 

W(R) = (V(r)) = (V(R + a)) = ^1 J d 3 s e~& V(R + s) . (27) 



The standard deviation a is given by Equation ()23|) , Note that the effective potential W 
depends on K, since a depends on K. 

To summarize, we propose the approximation 

Z « C 3 J d 3 a |lIC„ 6 / d 6 a n ) exp , (28) 

where 

^ = ^C n 2 |a n | 2 + -y o drW(iJ(r)). (29) 
The most drastic approximation is of course to take K = 0, see [S]. Then we get 



« = M, (30) 
V 12m 



and 



Z ra ( T7T^72 ) / d "o exp(-^W(oo)) ■ (31) 



/ m 
\2irf3h 2 

The present version of the partial averaging method is simpler than the original one |10[ 
I1H 1121 ITS] , in that the standard deviation a is taken to be constant. The method has been 
used previously for computing the propagator, and not directly the partition function. Then 
a has to vary along the path, since it must vanish at the end points. 

4 Example 1: The harmonic oscillator 

If V is a harmonic oscillator potential, 

V(r) = - mwV , (32) 

then the averaged potential W is just V plus a constant, 

W(r) = V(r) + - muj 2 a 2 . (33) 

The addition to the potential contributes a multiplicative factor in the partition function, and 
the resulting approximation is 



4tt 2 \ 6 

Another way to obtain the same approximation is to set 

for n > X, in Equation (|87|) with i? = 0, this is valid when we choose K large enough that 

« 1 • (38) 
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5 Example 2: The Coulomb potential 

Consider now two particles of masses mi, 7713 and charges qi,(}2, interacting by the Coulomb 
potential 

V = V( ri ,r 2 ) = - ^ (37) 

47reo — T2I 

We write Fourier expansions for the paths of both particles, 

00 00 
r i( T ) = X! a i» ei ^ I > r 2 (r)= X a 2» ei ^' ( 38 ) 

n=— 00 ra=— 00 

The partition function is 



MfP)7 d Wnp^) /"-H- 1 ' (39) 



where 



S = J dr f-m 1 (f 1 (r)) 2 + -m 2 (r'2(r)) 2 + y(ri(T),r2(r))J 

47T 2 °° fP h 

= X n2 ( m i l a i«! 2 + m 2 l a 2n| 2 ) + / dr V(t(t)) . (40) 
The potential depends only on the relative position 

r(r) =ri (r)-r 2 (r)= X fl ne , (41) 



where a n = a\ n — a 2n . The real and imaginary parts of the x, y, z components of the Fourier 
coefficients a„ have mean values zero and standard deviations 



= (|a lra | 2 ) + (|a 2n | 2 ) = _h_ J_ 
1 6 2nvr V 2m ' 



where m is the reduced mass, 



111 

- = — + —. 43 

m mi m 2 

Like in the one particle case, we define R(t) by an equation of the same form as Equation (|19|) . 
We integrate explicitly over the Fourier coefficients a\ n and a 2n up to n = K, and we do the 
remaining integrations approximately by averaging the potential as in Equation l|27|). The 
averaged Coulomb potential is 

W(r) = erff . (44) 

It equals the Coulomb potential in the limit r — > oo, but is nonsingular at the origin. The 
standard deviation cr is defined as in Equation (|2*3*|) . now with m as the reduced mass. The 
effect of the averaging is a multiplication by the error function, defined as 

erf(x) = -= / X du e~ u * . (45) 



vr Jo 
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K 


1 


2 


3 


4 


5 


o 


0.211 










1 


0.265 


1.068 








2 


0.273 


1.110 


2.032 






3 


0.276 


1.123 


2.059 


3.018 




4 


0.277 


1.128 


2.070 


3.036 


4.011 


5 


0.278 


1.131 


2.075 


3.045 


4.024 


10 


0.278 


1.135 


2.083 


3.058 


4.042 


20 


0.279 


1.136 


2.086 


3.061 


4.047 


50 


0.279 


1.137 


2.086 


3.062 


4.048 



Table 1: Zeros of the correlation functions fx(u), with < u < 1/2, multiplied by 2K + 1. 

6 Numerical computation 

A numerical estimate of the action Sa with the averaged potential W is the completely 
discretized action Sd defined by 

a K a 2K 



E^KP + ^yE^), (46) 



n=l 



where the positions rj are given by the Fourier coefficients according to Equation (|llj) . 

There are at least three arguments in favour of choosing exactly 2K + 1 evaluation points 
for the action integral of the potential. One is that this replacement of the integral by a sum is 
exact for a constant, linear or quadratic potential. Another argument is that the real Fourier 
coefficient ao and the K complex Fourier coefficients oi, . . . , are just what is needed to 
fix the 2K + 1 positions rj . 

The third argument is less obvious. In fact, our justification of the averaging procedure 
defining W suggests that the evaluation points Tj should be chosen in such a way that the 
covariances (s(rj) ■ s(rjt)), given by Equation (|21|) . are small. These covariances are propor- 
tional to the function fx(u) given in Equation Q22|) . which is symmetric about u = 1/2, and 
has K + 1 zeros between and 1/2. Of these zeros, K are close to the values Uj = j/(2K + 1), 
for j = 1,2, ... ,K, as Tabled shows. Thus, with 2K + 1 equally spaced points Tj we have 
that (s(Tj) • s(rfc)) k, for j ^ k. 

This does not necessarily forbid us to use for example twice as many evaluation points. 
We may define At = ph/(2K + 1) and 

2K 2K 

S a = Ar^^WAr)) , S b = AT^WiRdj + i)Ar)) . (47) 

3=0 j=0 

S a _ N Sg+Sfr 

But if we do so, we should perhaps compute (e s + e n )/2 rather than e za . 

To summarize again, the numerical approximation proposed here is based on evaluation 
of the following integral, where Sd is the discrete action defined in Equation (|46|) . 

Z « C 3 J d 3 a (j[C*J d 6 a„,^ exp (" ^) • (48) 
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From the partition function we compute the expectation value of the energy as 

E = (H) = --^lnZ((3). (49) 

In the above approximate partition function there is (3 dependence in the coefficients C n and 
C n , and also in the averaged potential W. We get that 



Jd 3 a (n^i/dV) 


ex P (-^) 


I 
n 


dS D 

a/3 


/d 3 a 




any 


exp(— 


n j 





E « v ' y + ^ y r ; 7 (50) 



i 3s D i #^ 2| l2 i ^ / T , . „ a 



and 

Note that Equation (|5()j) may be rewritten as 

Here the Fourier coefficients are present only in the kinetic part of the discrete action Sr>. 
In general we have that 

e%i w u-e%i; w u-ii; w u- (53) 

Hence we get, in the example of the harmonic oscillator potential, Equation (j.33j) . 

^W(r) = -mwV. (54) 
And in the example of the Coulomb potential, Equation (|44j). 

/3 £ W(r) = exp ( ^4) . (55) 

5/5 V ' 4vreo(^^) \ 2cJ / 

In the approximate expression for (H) = (T) + (V) it is not immediately obvious which 
contributions represent kinetic and potential energy, respectively. In order to identify the 
different terms, we should define 

Z(p 1 ,p 2 ) = Tre-^ T ~^ v , (56) 

and use that, e.g., 

(V) = --£-]nZ(p 1 ,p 2 ) . (57) 

This formula holds because 



9 e -/3iT-/3 2 V _ _ f 1 dA e -A(/3iT+/3 2 V) ^-(i-AJOSiT+jSaV) ^ ( 5 g) 



2 
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and hence 




Tr(Ve-^ T ~^ v ) . 



(59) 



The somewhat surprising conclusion is that 



(V) 



(ngo/dV 3 -)exp(-^) 



(60) 



(ngo/dV,) exp 




whereas (T) is all the rest of the right hand side of Equation ((50(1 or Equation 1152(1 . 

For bound states of any number of particles, with the Coulomb interaction, the virial 
theorem states that 2(T) + (V) = 0. It gives a good check on numerical results for bound 
states, if one computes both (T) and (V). It may also be used to (potentially) improve the 
precision of computed energies, since it implies for example that 



The statistical error with which the two expectation values (T) and (V) are computed in a 
Monte Carlo simulation will in general not be the same, hence one may use whichever value 
has the smallest error. 

7 Regularization 

For simplicity, we have so far neglected the fact that the partition function is not mathemat- 
ically well defined for a system in an infinite volume when, for example, the potential goes to 
zero at infinity, like the Coulomb potential. In our present context, the problem is that the 
integral over ao diverges. In practice, when the integral is computed by some Monte Carlo 
method using a random walk algorithm of the Metropolis type, the divergence means that 
there is a finite probability of walking away to infinity, where the potential vanishes. This 
may be no problem in practice, because the divergence may be so improbable that it will 
never happen in the Monte Carlo simulation. Nevertheless, one may like to introduce some 
kind of regularization which makes the partition function well defined. 

A convenient regularization method in our case is to add to the Hamiltonian an extra 
harmonic oscillator potential 



with a suitably chosen angular frequency ujq. The Fourier expansion of Equation (J5J) implies 
that 



E = (T) + (V} = -(T} = Y- 



(61) 




(62) 




(63) 



Hence, Equation Q is modified to read 




(64) 
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where 



A? = ^ , ^ = + /W = 4( - 2 t:> 2 - , (65) 



for n = 1, 2, . . ., and 

V = . 

2vr 

We define also D n = D n / ^/tt for n = 0, 1, 2, Since 



(66) 



n oo ^ 2 i oo o i 

°0 TT u n _ 1 TT n _ 1 

Do D 2 " /3^ ^ n 2 + v 2 ~ 2 sinh^vr) ' 1 ' 



we may rewrite Equation © as 

Z = Z 5 3 1 d 3 a f] ( 5 n / d V ) exp 

n=l 



(68) 



Z = Z (f3) = - . i3 ~ . (69) 



where Zo is the partition function of the three dimensional harmonic oscillator with angular 
frequency u>o, 

1 

! sinh 3 (z/7r) 

This expression for the partition function Z = Z{(5) is mathematically well defined, when 
the harmonic oscillator potential V$ is included in addition to the potential V, so that Equa- 
tion ((OH) holds. 

A natural way to interpret Equation (|68[) is that the Fourier coefficients a n are Gaussian 
random variables with mean zero and standard deviations a n = l/(y/2 D n ). Our derivation of 
how to replace the potential V by an averaged potential W, goes through with little change. 
The most important change is that the denominator n 2 in Equation 1)21)1 has to be replaced 
by n 2 + v 2 , and hence the correlation function fx is replaced by a function fx „ which is still 
periodic with period 1, 

, . ^ cos(2n7ru) 7r cosh(^7r(2n — 1)) 1 ^ cos(2n7rn) , , 

JkA u ) = 2^ — ~, — ~ = o — r T7 — ^ ~ — ~, — T~ ' v 70 ) 

n=K+l v ' n=l 

the last formula being valid for < u < 1. The zeros of /k^ between and 1/2 are even 
closer to the values u& = k/(2K + 1), for k = 1,2, . . . , K , than those of fx, as Table |2 shows, 
for the arbitrarily chosen value v = 10. 

The numerical computation now involves the following modified version of Equation 1)48)1 , 

Z « Z A? | d 3 a (il 5 n / d 6 a„) ex p("x) ' 
where we use also a modified definition of the discrete action Sd, 

71=0 J=0 

if a 2K" 

= £ Q 2 |a n | 2 + J2 (Vo^) + W( rj )) . (72) 

n=l j=0 
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K 


1 


2 


3 


4 


5 


o 


0.055 










1 


0.113 


1.003 








2 


0.151 


1.011 


2.002 






3 


0.177 


1.022 


2.008 


3.002 




4 


0.197 


1.035 


2.015 


3.007 


4.002 


5 


0.211 


1.048 


2.021 


3.012 


4.006 


10 


0.249 


1.093 


2.050 


3.033 


4.023 


20 


0.269 


1.122 


2.073 


3.051 


4.038 


50 


0.277 


1.134 


2.084 


3.060 


4.046 



Table 2: Zeros of the correlation functions //^(n), with < u < 1/2, multiplied by 2K + 1. 
The table is for v = 10. 

The standard deviation a to be used in the definition of the averaged potential W, will now 
be given by the formula 

° ~ 2^m lKA ' ~ 2^m { 2v 2v* ^ n 2 + ) ' { ' 6) 

Using Equation (|71|). we compute the total energy (E) = (T) + (V) + (Vo) as 

a jgg+i) 6^ (n, 2 io/dS)ex P (-^)^ 

< £ > - '» z — sr- + t /i '"' (0) + (S|57^M=¥) ' (74) 

Compare this to Equation (|52|). The modified version of Equation (|51(1 is the following, 



„ fl/? ^ c ^ 2 + wri^X° (ri)+w{ri)+P ^ w{ri) l- (75) 

Equation (|53|) gets modified as follows, 

P iL W(r) = (l + u-?-\n f K>v (0)) °- #- W(r) . (76) 



0(3 w V Si/ ' ' v 7 2 0cj 
In order to calculate separately the expectation values (T), (V), and (Vo), we should define 

Z = Z(/3 ,/3i,/3 2 ) = Tre-^ T -^ v -P oVo , (77) 

and keep track of the three parameters Po,{h>(h before setting them all equal to (3. This 
gives for (V) a formula exactly like Equation IjBOJI . An easier way to compute (Vo) is to note 
that 

lv ._ o,„ 9 ]nZ _3^ (Tig, J Ah) lm 

< F °> - " 25 as 1,1 z ~ T /t>(0) + (nftjasHf ^ ' ( ' 
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Here we have that 

w o dS D 1 ^ / . . u; 3 A 

2^ = (2FTT) g r (rj) Ta^ ' (79) 

with 



Once we know (E), (V) and {Vq}, we know also (T) = (E) — (V) — {Vq). The identification 
of the various contributions to the total energy (E) is seen to be not entirely trivial. 
If V is taken to be the Coulomb potential, the virial theorem gives now that 

2(T) + (V) - 2{V ) = . (81) 

This holds for any number of particles. It provides a check on numerical results, and it may 
be used to compute the total energy, including the regulator potential, as 

E = (T) + (V) + (V ) = ^ + 2(V ) . (82) 
The energy including the Coulomb potential but excluding the regulator potential, is 

E C = (T) + (V} = ^± + (V ) . (83) 



8 Numerical test results 

Table [3 presents numerical results for the ground state energy of the hydrogen atom, for 
comparison with the exact value of —13.598 eV. All results are for a temperature of 15 000 K. 

The Monte Carlo method was used with a standard Metropolis algorithm. In each Monte 
Carlo step, one point to be updated is chosen randomly among the 2K + 1 points on the 
discrete path, then a random step is generated and either accepted or rejected depending on 
the change in the discrete action, ASA. If ASd < 0, the step is accepted. If ASd > — hlau, 
with u a uniform random variable between and 1, the step is rejected. The optimization 
of the Monte Carlo strategy was not considered, but is of course an important problem. In 
fact, the naive approach of updating one point at a time has a disastrously slow convergence 
when more than about one hundred time steps are used. 

The main computational cost of updating one point is computing the change in the Fourier 
components, this takes approximately 2K+1 floating point operations. If more than one point 
is updated in each step, one may choose e.g. 2K + 1 = 3 n for some power n, and then use the 
fast Fourier transform with base 3. 

By far the highest statistics, 10 10 MC steps, was run for the entry with 201 time steps. In 
this case, statistical uncertainties are given in the table, and the values found for the ground 
state energy are consistent with the exact value, within the uncertainties of less than one 
per cent. A number of time steps of the order of 50 may give sufficient accuracy for many 
purposes. 

It is noteworthy that the statistical error in the direct estimate of the energy, (T) + (V), 
is half the separate errors in (T) and (V). In fact, (T) + (V) is seen to be systematically 
closer to the exact energy than the estimate (Vo) + (V)/2 obtained from the virial theorem 
by elimination of the kinetic energy. 
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No. of points 


Ener 


eies 


Regularization 


Virial 


Coulomb energy 


2K + 1 


(T) 


(V) 




(V a ) 
\ r u/ 


A 


(T) + (V) 


(V n ) + (V)/2 

\* U/ I \ v / / ^ 


21 


12.157 


-25.486 


o 


0.000 


0.014 


-12.729 


-12.743 


41 


13.063 


-26 098 


o 


0.000 


0.014 


-13.035 


-13.049 


1 


6.538 


-13.193 


1 


0.146 


-0.205 


-6.655 


-6.450 


11 


11.676 


-23 680 


1 


0.064 


-0.228 


-12.004 


-11.776 


21 


12.513 


-25.351 


1 


0.061 


-0.223 


-12.838 


-12.615 


41 


13.134 


-26.747 


1 


0.054 


-0.293 


-13.613 


-13.320 


101 


13.890 


-27.613 


1 


0.050 


0.034 


-13.723 


-13.757 


201 


13.489 


-26.991 


1 


0.056 


-0.063 


-13.502 


-13.439 




±0.141 


±0.181 




±0.002 


±0.064 


±0.070 


±0.093 


1 


13.417 


-20.191 


10 


3.456 


-0.135 


-6.774 


-6.640 


11 


17.728 


-28.538 


10 


3.599 


-0.140 


-10.810 


-10.670 


21 


18.520 


-30.363 


10 


3.461 


-0.122 


-11.843 


-11.721 


41 


18.958 


-31.456 


10 


3.381 


-0.151 


-12.498 


-12.347 



Table 3: Estimates of the hydrogen ground state energy. All values tabulated are in eV. The 
quantity A = (T) — (Vq) + {V)/2 should be zero, by the virial theorem. The last two columns 
should be compared to the exact value of —13.598 eV. See comments in the text. 
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A Example: Quadratic Lagrangian 

In the case of a particle of electric charge q moving in an electromagnetic vector potential 
A(r), the imaginary time action has also an imaginary part, 

S = dr Q m (r(r)) 2 + V{r{r)) + i q f(r) • A(r(r))j . (84) 

Note that the contribution from the vector potential is gauge invariant, because we integrate 
over a closed path. The partition function can be computed exactly by the path integral for 
example when we have an isotropic harmonic oscillator external potential of angular frequency 
(jj, and a magnetic field of constant flux density B, so that 

V{r)= l -mu 2 r\ A(r) = . (85) 

Then 

| = /W |a | 2 + £ ((C n 2 + /W) K| 2 + B ■ (a n x <)) , (86) 
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and the partition function is, with B = \B\, 



Z 



n 



1 + 



2nir) 



1 + 



V 2nvr 



)T+( 



^ iggj 

2nirm 



(87) 



The energy spectrum is of course well known. An energy eigenvalue is given by quantum 
numbers j, k, £ = 0, 1, 2, . . . as 



where 



Hence, 



E 



■j,k,e 



l 



l 



l 



Iuj 2 + 



f\qB\Y ± \qB\ 



V 2 



2m 



1 



j,M 8 sinh ( ^ sinh ^ ^ sinh^ 

Equation (fHTj) gives a product representation of this function. 



(90) 
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